/* sstebz.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
 */

#include "Lapack.h"
#include <cmath>
#include "../Profiler.h"

using namespace std;

extern "C" {
	int lsame_(char *, char *);
	int xerbla_(char *, int *);
}

extern int ilaenv_(int *, const char *, const char *, int *, int *, int *, int *);
extern double slamch_(char *);

/* Table of constant values */

static int c__1 = 1;
static int c_n1 = -1;
static int c__3 = 3;
static int c__2 = 2;
static int c__0 = 0;

enum Range {
	RANGE_NONE = 0,
	RANGE_ALL,
	RANGE_VALUE,
	RANGE_INDEX,
};

enum Order {
	ORDER_NONE = 0,
	ORDER_ENTIRE_MATRIX,
	ORDER_BY_BLOCK,
};

struct Interval {
	float low;
	float up;
};

int Lapack::sstebz(char *range, char *order, int *n, float *vl,
		float *vu, int *il, int *iu, float *abstol, float *d, float *e,
		int *m, int *nsplit, float *w, int *iblock, int *
		isplit, float *work, int *iwork, int *info) 
{
	Measure measure("sstebz");
	
	/* Local variables */
	int j, ib, jb, ie, je, nb;
	float gl;
	int im, in;
	float gu;
	int iw;
	float wl, wu;
	int nwl;
	float ulp, wlu, wul;
	int nwu;
	float tmp1, tmp2;
	int iend, ioff, iout, itmp1, jdisc;
	int iinfo;
	float atoli;
	int iwoff;
	float bnorm;
	int itmax;
	float wkill, rtoli, tnorm;
	int ibegin, irange, idiscl;
	float safemn;
	int idumma[1];
	int idiscu;
	int iorder;
	int ncnvrg;
	float pivmin;
	int toofew;


	/*  -- LAPACK routine (version 3.2) -- */
	/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
	/*     November 2006 */
	/*     8-18-00:  Increase FUDGE factor for T3E (eca) */

	/*     .. Scalar Arguments .. */
	/*     .. */
	/*     .. Array Arguments .. */
	/*     .. */

	/*  Purpose */
	/*  ======= */

	/*  SSTEBZ computes the eigenvalues of a symmetric tridiagonal */
	/*  matrix T.  The user may ask for all eigenvalues, all eigenvalues */
	/*  in the half-open interval (VL, VU], or the IL-th through IU-th */
	/*  eigenvalues. */

	/*  To avoid overflow, the matrix must be scaled so that its */
	/*  largest element is no greater than overflow**(1/2) * */
	/*  underflow**(1/4) in absolute value, and for greatest */
	/*  accuracy, it should not be much smaller than that. */

	/*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
	/*  Matrix", Report CS41, Computer Science Dept., Stanford */
	/*  University, July 21, 1966. */

	/*  Arguments */
	/*  ========= */

	/*  RANGE   (input) CHARACTER*1 */
	/*          = 'A': ("All")   all eigenvalues will be found. */
	/*          = 'V': ("Value") all eigenvalues in the half-open interval */
	/*                           (VL, VU] will be found. */
	/*          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */
	/*                           entire matrix) will be found. */

	/*  ORDER   (input) CHARACTER*1 */
	/*          = 'B': ("By Block") the eigenvalues will be grouped by */
	/*                              split-off block (see IBLOCK, ISPLIT) and */
	/*                              ordered from smallest to largest within */
	/*                              the block. */
	/*          = 'E': ("Entire matrix") */
	/*                              the eigenvalues for the entire matrix */
	/*                              will be ordered from smallest to */
	/*                              largest. */

	/*  N       (input) INTEGER */
	/*          The order of the tridiagonal matrix T.  N >= 0. */

	/*  VL      (input) REAL */
	/*  VU      (input) REAL */
	/*          If RANGE='V', the lower and upper bounds of the interval to */
	/*          be searched for eigenvalues.  Eigenvalues less than or equal */
	/*          to VL, or greater than VU, will not be returned.  VL < VU. */
	/*          Not referenced if RANGE = 'A' or 'I'. */

	/*  IL      (input) INTEGER */
	/*  IU      (input) INTEGER */
	/*          If RANGE='I', the indices (in ascending order) of the */
	/*          smallest and largest eigenvalues to be returned. */
	/*          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
	/*          Not referenced if RANGE = 'A' or 'V'. */

	/*  ABSTOL  (input) REAL */
	/*          The absolute tolerance for the eigenvalues.  An eigenvalue */
	/*          (or cluster) is considered to be located if it has been */
	/*          determined to lie in an interval whose width is ABSTOL or */
	/*          less.  If ABSTOL is less than or equal to zero, then ULP*|T| */
	/*          will be used, where |T| means the 1-norm of T. */

	/*          Eigenvalues will be computed most accurately when ABSTOL is */
	/*          set to twice the underflow threshold 2*SLAMCH('S'), not zero. */

	/*  D       (input) REAL array, dimension (N) */
	/*          The n diagonal elements of the tridiagonal matrix T. */

	/*  E       (input) REAL array, dimension (N-1) */
	/*          The (n-1) off-diagonal elements of the tridiagonal matrix T. */

	/*  M       (output) INTEGER */
	/*          The actual number of eigenvalues found. 0 <= M <= N. */
	/*          (See also the description of INFO=2,3.) */

	/*  NSPLIT  (output) INTEGER */
	/*          The number of diagonal blocks in the matrix T. */
	/*          1 <= NSPLIT <= N. */

	/*  W       (output) REAL array, dimension (N) */
	/*          On exit, the first M elements of W will contain the */
	/*          eigenvalues.  (SSTEBZ may use the remaining N-M elements as */
	/*          workspace.) */

	/*  IBLOCK  (output) INTEGER array, dimension (N) */
	/*          At each row/column j where E(j) is zero or small, the */
	/*          matrix T is considered to split into a block diagonal */
	/*          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which */
	/*          block (from 1 to the number of blocks) the eigenvalue W(i) */
	/*          belongs.  (SSTEBZ may use the remaining N-M elements as */
	/*          workspace.) */

	/*  ISPLIT  (output) INTEGER array, dimension (N) */
	/*          The splitting points, at which T breaks up into submatrices. */
	/*          The first submatrix consists of rows/columns 1 to ISPLIT(1), */
	/*          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
	/*          etc., and the NSPLIT-th consists of rows/columns */
	/*          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
	/*          (Only the first NSPLIT elements will actually be used, but */
	/*          since the user cannot know a priori what value NSPLIT will */
	/*          have, N words must be reserved for ISPLIT.) */

	/*  WORK    (workspace) REAL array, dimension (4*N) */

	/*  IWORK   (workspace) INTEGER array, dimension (3*N) */

	/*  INFO    (output) INTEGER */
	/*          = 0:  successful exit */
	/*          < 0:  if INFO = -i, the i-th argument had an illegal value */
	/*          > 0:  some or all of the eigenvalues failed to converge or */
	/*                were not computed: */
	/*                =1 or 3: Bisection failed to converge for some */
	/*                        eigenvalues; these eigenvalues are flagged by a */
	/*                        negative block number.  The effect is that the */
	/*                        eigenvalues may not be as accurate as the */
	/*                        absolute and relative tolerances.  This is */
	/*                        generally caused by unexpectedly inaccurate */
	/*                        arithmetic. */
	/*                =2 or 3: RANGE='I' only: Not all of the eigenvalues */
	/*                        IL:IU were found. */
	/*                        Effect: M < IU+1-IL */
	/*                        Cause:  non-monotonic arithmetic, causing the */
	/*                                Sturm sequence to be non-monotonic. */
	/*                        Cure:   recalculate, using RANGE='A', and pick */
	/*                                out eigenvalues IL:IU.  In some cases, */
	/*                                increasing the PARAMETER "FUDGE" may */
	/*                                make things work. */
	/*                = 4:    RANGE='I', and the Gershgorin interval */
	/*                        initially used was too small.  No eigenvalues */
	/*                        were computed. */
	/*                        Probable cause: your machine has sloppy */
	/*                                        floating-point arithmetic. */
	/*                        Cure: Increase the PARAMETER "FUDGE", */
	/*                              recompile, and try again. */

	/*  Internal Parameters */
	/*  =================== */

	/*  RELFAC  REAL, default = 2.0e0 */
	/*          The relative tolerance.  An interval (a,b] lies within */
	/*          "relative tolerance" if  b-a < RELFAC*ulp*max(|a|,|b|), */
	/*          where "ulp" is the machine precision (distance from 1 to */
	/*          the next larger floating point number.) */

	/*  FUDGE   REAL, default = 2 */
	/*          A "fudge factor" to widen the Gershgorin intervals.  Ideally, */
	/*          a value of 1 should work, but on machines with sloppy */
	/*          arithmetic, this needs to be larger.  The default for */
	/*          publicly released versions should be large enough to handle */
	/*          the worst machine around.  Note that this has no effect */
	/*          on accuracy of the solution. */

	/*  ===================================================================== */

	/*     .. Parameters .. */
	/*     .. */
	/*     .. Local Scalars .. */
	/*     .. */
	/*     .. Local Arrays .. */
	/*     .. */
	/*     .. External Functions .. */
	/*     .. */
	/*     .. External Subroutines .. */
	/*     .. */
	/*     .. Intrinsic Functions .. */
	/*     .. */
	/*     .. Executable Statements .. */

	/* Parameter adjustments */
	--iwork;
	--work;
	--isplit;
	--iblock;
	--w;
	--e;
	--d;

	/* Function Body */
	*info = 0;

	/*     Decode RANGE */

	if (lsame_(range, "A")) {
		irange = RANGE_ALL;
	} else if (lsame_(range, "V")) {
		irange = RANGE_VALUE;
	} else if (lsame_(range, "I")) {
		irange = RANGE_INDEX;
	} else {
		irange = RANGE_NONE;
	}

	/*     Decode ORDER */

	if (lsame_(order, "B")) {
		iorder = ORDER_BY_BLOCK;
	} else if (lsame_(order, "E")) {
		iorder = ORDER_ENTIRE_MATRIX;
	} else {
		iorder = ORDER_NONE;
	}

	/*     Check for Errors */

	if (irange <= 0) {
		*info = -1;
	} else if (iorder <= 0) {
		*info = -2;
	} else if (*n < 0) {
		*info = -3;
	} else if (irange == RANGE_VALUE) {
		if (*vl >= *vu) {
			*info = -5;
		}
	} else if (irange == RANGE_INDEX && (*il < 1 || *il > max(1, *n))) {
		*info = -6;
	} else if (irange == RANGE_INDEX && (*iu < min(*n, *il) || *iu > *n)) {
		*info = -7;
	}

	if (*info != 0) {
		int temp = -(*info);
		xerbla_("SSTEBZ", &temp);
		return 0;
	}

	/*     Initialize error flags */

	*info = 0;
	ncnvrg = false;
	toofew = false;

	/*     Quick return if possible */

	*m = 0;
	if (*n == 0) {
		return 0;
	}

	/*     Simplifications: */

	if (irange == RANGE_INDEX && *il == 1 && *iu == *n) {
		irange = RANGE_ALL;
	}

	/*     Get machine constants */
	/*     NB is the minimum vector length for vector bisection, or 0 */
	/*     if only scalar is to be done. */

	safemn = slamch_("S");
	ulp = slamch_("P");
	rtoli = ulp * 2.f;
	nb = ilaenv_(&c__1, "SSTEBZ", " ", n, &c_n1, &c_n1, &c_n1);
	if (nb <= 1) {
		nb = 0;
	}

	/*     Special Case when N=1 */

	if (*n == 1) {
		*nsplit = 1;
		isplit[1] = 1;
		if (irange == RANGE_VALUE && (*vl >= d[1] || *vu < d[1])) {
			*m = 0;
		} else {
			w[1] = d[1];
			iblock[1] = 1;
			*m = 1;
		}
		return 0;
	}

	/*     Compute Splitting Points */

	*nsplit = 1;
	work[*n] = 0.f;
	pivmin = 1.f;

	for (j = 2; j <= *n; ++j) {
		tmp1 = e[j - 1] * e[j - 1];
		if (abs(d[j] * d[j - 1]) * (ulp * ulp) + safemn > tmp1) {
			isplit[*nsplit] = j - 1;
			++(*nsplit);
			work[j - 1] = 0.f;
		} else {
			work[j - 1] = tmp1;
			pivmin = max(pivmin, tmp1);
		}
	}
	isplit[*nsplit] = *n;
	pivmin *= safemn;

	/*     Compute Interval and ATOLI */

	if (irange == RANGE_INDEX) {

		/*        RANGE='I': Compute the interval containing eigenvalues */
		/*                   IL through IU. */

		/*        Compute Gershgorin interval for entire (split) matrix */
		/*        and use it as the initial interval */

		gu = d[1];
		gl = d[1];
		tmp1 = 0.f;

		for (j = 1; j <= *n - 1; ++j) {
			tmp2 = sqrt(work[j]);
			gu = max(gu, d[j] + tmp1 + tmp2);
			gl = min(gl, d[j] - tmp1 - tmp2);
			tmp1 = tmp2;
		}

		gu = max(gu, d[*n] + tmp1);
		gl = min(gl, d[*n] - tmp1);
		tnorm = max(abs(gl), abs(gu));
		gl = gl - tnorm * 2.1f * ulp * *n - pivmin * 4.2000000000000002f;
		gu = gu + tnorm * 2.1f * ulp * *n + pivmin * 2.1f;

		/*        Compute Iteration parameters */

		itmax = (int) ((log(tnorm + pivmin) - log(pivmin)) / log(2.f)) + 2;
		if (*abstol <= 0.f) {
			atoli = ulp * tnorm;
		} else {
			atoli = *abstol;
		}

		work[*n + 1] = gl;
		work[*n + 2] = gl;
		work[*n + 3] = gu;
		work[*n + 4] = gu;
		work[*n + 5] = gl;
		work[*n + 6] = gu;
		iwork[1] = -1;
		iwork[2] = -1;
		iwork[3] = *n + 1;
		iwork[4] = *n + 1;
		iwork[5] = *il - 1;
		iwork[6] = *iu;

		slaebz(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, &pivmin,
				&d[1], &e[1], &work[1], &iwork[5], &work[*n + 1], &work[*n + 5],
				&iout, &iwork[1], &w[1], &iblock[1], &iinfo);

		if (iwork[6] == *iu) {
			wl = work[*n + 1];
			wlu = work[*n + 3];
			nwl = iwork[1];
			wu = work[*n + 4];
			wul = work[*n + 2];
			nwu = iwork[4];
		} else {
			wl = work[*n + 2];
			wlu = work[*n + 4];
			nwl = iwork[2];
			wu = work[*n + 3];
			wul = work[*n + 1];
			nwu = iwork[3];
		}

		if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
			*info = 4;
			return 0;
		}
	} else {

		/*        RANGE='A' or 'V' -- Set ATOLI */
		tnorm = max(abs(d[1]) + abs(e[1]), abs(d[*n]) + abs(e[*n - 1]));
		for (j = 2; j <= *n - 1; ++j) {
			tnorm = max(tnorm, abs(d[j]) + abs(e[j - 1]) + abs(e[j]));
		}

		if (*abstol <= 0.f) {
			atoli = ulp * tnorm;
		} else {
			atoli = *abstol;
		}

		if (irange == RANGE_VALUE) {
			wl = *vl;
			wu = *vu;
		} else {
			wl = 0.f;
			wu = 0.f;
		}
	}

	/*     Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. */
	/*     NWL accumulates the number of eigenvalues .le. WL, */
	/*     NWU accumulates the number of eigenvalues .le. WU */

	*m = 0;
	iend = 0;
	*info = 0;
	nwl = 0;
	nwu = 0;

	for (jb = 1; jb <= *nsplit; ++jb) {
		ioff = iend;
		ibegin = ioff + 1;
		iend = isplit[jb];
		in = iend - ioff;

		if (in == 1) {

			/*           Special Case -- IN=1 */

			if (irange == RANGE_ALL || wl >= d[ibegin] - pivmin) {
				++nwl;
			}
			if (irange == RANGE_ALL || wu >= d[ibegin] - pivmin) {
				++nwu;
			}
			if (irange == RANGE_ALL || wl < d[ibegin] - pivmin && wu >= d[ibegin]
					- pivmin) {
				++(*m);
				w[*m] = d[ibegin];
				iblock[*m] = jb;
			}
		} else {

			/*           General Case -- IN > 1 */

			/*           Compute Gershgorin Interval */
			/*           and use it as the initial interval */

			gu = d[ibegin];
			gl = d[ibegin];
			tmp1 = 0.f;

			for (j = ibegin; j <= iend - 1; ++j) {
				tmp2 = abs(e[j]);
				gu = max(gu, d[j] + tmp1 + tmp2);
				gl = min(gl, d[j] - tmp1 - tmp2);
				tmp1 = tmp2;
			}

			gu = max(gu, d[iend] + tmp1);
			gl = min(gl, d[iend] - tmp1);
			bnorm = max(abs(gl), abs(gu));
			gl = gl - bnorm * 2.1f * ulp * in - pivmin * 2.1f;
			gu = gu + bnorm * 2.1f * ulp * in + pivmin * 2.1f;

			/*           Compute ATOLI for the current submatrix */

			if (*abstol <= 0.f) {
				atoli = ulp * max(abs(gl), abs(gu));
			} else {
				atoli = *abstol;
			}

			if (irange > 1) {
				if (gu < wl) {
					nwl += in;
					nwu += in;
					goto L70;
				}
				gl = max(gl, wl);
				gu = min(gu, wu);
				if (gl >= gu) {
					goto L70;
				}
			}

			/*           Set Up Initial Interval */

			work[*n + 1] = gl;
			work[*n + in + 1] = gu;
			slaebz(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, &
					pivmin, &d[ibegin], &e[ibegin], &work[ibegin], idumma, &
					work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
					w[*m + 1], &iblock[*m + 1], &iinfo);

			nwl += iwork[1];
			nwu += iwork[in + 1];
			iwoff = *m - iwork[1];

			/*           Compute Eigenvalues */

			itmax = (int) ((log(gu - gl + pivmin) - log(pivmin)) / log(
					2.f)) + 2;
			slaebz(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, &
					pivmin, &d[ibegin], &e[ibegin], &work[ibegin], idumma, &
					work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
					&w[*m + 1], &iblock[*m + 1], &iinfo);

			/*           Copy Eigenvalues Into W and IBLOCK */
			/*           Use -JB for block number for unconverged eigenvalues. */

			for (j = 1; j <= iout; ++j) {
				tmp1 = (work[j + *n] + work[j + in + *n]) * .5f;

				/*              Flag non-convergence. */

				if (j > iout - iinfo) {
					ncnvrg = true;
					ib = -jb;
				} else {
					ib = jb;
				}
				int count = iwork[j + in] + iwoff;
				for (je = iwork[j] + 1 + iwoff; je <= count; ++je) {
					w[je] = tmp1;
					iblock[je] = ib;
				}
			}

			*m += im;
		}
L70:
		;
	}

	/*     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */
	/*     If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */

	if (irange == RANGE_INDEX) {
		im = 0;
		idiscl = *il - 1 - nwl;
		idiscu = nwu - *iu;

		if (idiscl > 0 || idiscu > 0) {
			for (je = 1; je <= *m; ++je) {
				if (w[je] <= wlu && idiscl > 0) {
					--idiscl;
				} else if (w[je] >= wul && idiscu > 0) {
					--idiscu;
				} else {
					++im;
					w[im] = w[je];
					iblock[im] = iblock[je];
				}
			}
			*m = im;
		}
		if (idiscl > 0 || idiscu > 0) {

			/*           Code to deal with effects of bad arithmetic: */
			/*           Some low eigenvalues to be discarded are not in (WL,WLU], */
			/*           or high eigenvalues to be discarded are not in (WUL,WU] */
			/*           so just kill off the smallest IDISCL/largest IDISCU */
			/*           eigenvalues, by simply finding the smallest/largest */
			/*           eigenvalue(s). */

			/*           (If N(w) is monotone non-decreasing, this should never */
			/*               happen.) */

			if (idiscl > 0) {
				wkill = wu;
				for (jdisc = 1; jdisc <= idiscl; ++jdisc) {
					iw = 0;
					for (je = 1; je <= *m; ++je) {
						if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
							iw = je;
							wkill = w[je];
						}
					}
					iblock[iw] = 0;
				}
			}
			if (idiscu > 0) {

				wkill = wl;
				for (jdisc = 1; jdisc <= idiscu; ++jdisc) {
					iw = 0;
					for (je = 1; je <= *m; ++je) {
						if (iblock[je] != 0 && (w[je] > wkill || iw == 0)) {
							iw = je;
							wkill = w[je];
						}
					}
					iblock[iw] = 0;
				}
			}
			im = 0;
			for (je = 1; je <= *m; ++je) {
				if (iblock[je] != 0) {
					++im;
					w[im] = w[je];
					iblock[im] = iblock[je];
				}
			}
			*m = im;
		}
		if (idiscl < 0 || idiscu < 0) {
			toofew = true;
		}
	}

	/*     If ORDER='B', do nothing -- the eigenvalues are already sorted */
	/*        by block. */
	/*     If ORDER='E', sort the eigenvalues from smallest to largest */

	if (iorder == ORDER_ENTIRE_MATRIX && *nsplit > 1) {
		for (je = 1; je < *m; ++je) {
			ie = 0;
			tmp1 = w[je];
			for (j = je + 1; j <= *m; ++j) {
				if (w[j] < tmp1) {
					ie = j;
					tmp1 = w[j];
				}
			}

			if (ie != 0) {
				itmp1 = iblock[ie];
				w[ie] = w[je];
				iblock[ie] = iblock[je];
				w[je] = tmp1;
				iblock[je] = itmp1;
			}
		}
	}

	*info = 0;
	if (ncnvrg) {
		++(*info);
	}
	if (toofew) {
		*info += 2;
	}
	return 0;

	/*     End of SSTEBZ */

} /* sstebz_ */
